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We study the long time evolution and stationary speed distribution of N point particles in 2D 
moving under the action of an external field E, and undergoing elastic collisions with either a fixed 
f\j 1 periodic array of convex scatterers, or with virtual random scatterers. The total kinetic energy of 

the N particles is kept fixed by a Gaussian thermostat which induces an interaction between the 
O ■ particles. We show analytically and numerically that for weak fields this distribution is universal, 

' i.e. independent of the position or shape of the obstacles or the nature of the stochastic scattering. 

Our analysis is based on the existence of two time scales; the velocity directions become uniformized 
on: in times of order unity while the speeds change only on a time scale of 0(|E| ). 
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Our understanding of nonequilibrium stationary states (NESS) of multiparticle systems is very incomplete at 
present. In particular, there are no cases (we know of) where one has an explicit expression for the NESS of an 
interacting system of particles with positions and velocities. In this paper we derive an analytic expression for a 
non-trivial NESS having a certain universality. While our proof requires some technical assumptions the arguments 
are physically clear and convincing The results are furthermore checked by very extensive computer simulations. 

The model we consider is a variation of the Drude-Lorentz model of electrical conduction in two dimensions [lj]. The 
system consists of N particles (electrons) moving under the action of a constant external field E among a periodic 
array of fixed convex scatterers (Sinai billiard) with which they collide clastically. In order to produce a stationary 



■ current carrying state it is necessary to have a mechanism which will absorb the heat produced by the field E. This is 



modeled here by a Gaussian thermostat which keeps the total kinetic energy of the system constant [13[ . Thermostated 
systems are known to give results which are in accord with observations on real systems, c.f. the Gallavotti- Cohen 
£NJ . fluctuation relation 0, [ll| and the Ruelle results about heat conduction 14 1. 



The equations of motion for the system on the unit two dimensional torus, which corresponds to an infinite system 
with periodic scatterers and periodic initial conditions are, taking the mass of each particle to be one, 



I * = Vi i = l,...,N 

where 

N N 

j = 5>, t/ = $>*l 2 (2) 

i=l i=l 

and Ti is the "force" exerted on the zth particle by collisions with the fixed scatterers. These collisions only change 
the direction but not the speed of the particle. It is easy to see that due to the Gaussian thermostat, represented by 
the (E • J)vi/U term, the total kinetic energy of the system is constant, i.e., j^U = 0. 

This system was first introduced, for the case N = 1, by Moran and Hoover [bj ] where it was found numerically 
that the NESS had a fractal structure. This was shown rigorously in Q, where it was proven that this system has a 
unique singular SRB (Sinai-Ruelle-Bowen) measure for small E which satisfies Ohm's law. The N = 1 system was 
further investigated both numerically and analytically in Q Q ■ 

The multi-particle system N > 1 was first investigated numerically in These gave many tantalizing hints 
about the interesting nature of the NESS for this system resulting from the effective interaction between the particles 
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caused by the thermostat: if one particle increases (decreases) its speed due to the external field E the others have to 
decrease (increase) their speed to keep the total kinetic energy fixed. To study the NESS of this system analytically, 
a stochastic version of the dynamics was also introduced in |5| . In this stochastic model the collisions with the fixed 
obstacles are replaced by "virtual" collisions or scatterings. These collisions, like the fixed scatterers, conserve energy 
and tend to make the angular distribution uniform. It appeared numerically that the one particle marginal speed 
distribution of the stochastic and deterministic NESS for a symmetric billiard table (Table A in Figure [JJ were close 
for small E, but there was no argument of why this should be so. The accuracy of the simulations in Q was not very 
high, so it was assumed that these coincidences were approximate and valid only for the symmetric table. 

We recently revisited this problem using more computer power and new analytical methods. As a result we now 
have strong numerical and analytic evidence that the full NESS speed distribution is, in the limit E — > 0, exactly the 
same for the stochastic and deterministic models for all N. This implies ipso facto that in this limit the distribution is 
independent of the shape of the billiard table. It is a "universal" function whose exact shape we determine analytically 
below. Surprisingly the result seems to remain valid up to substantial values of |E|. Just how large E can be depends 
on the shape of the table; see Figures [TJ The new element in our analysis is the exploitation of a time scale separation 
which occurs for small |E|. 

Going beyond the limit E — > we find analytically the first order (in E) correction to the invariant distribution for 
a generalized version of the stochastic model. The expression we find can be extended to the deterministic billiard 
where, unlike the speed distribution, it depends on the shape of the table, but only through some properties which 
can be obtained from the N = 1 solution of the problem. We have checked this expression numerically by computing, 
with high precision, the invariant measure for Table B in Figure [TJ We also find analytically and verify numerically 
the asymptotic N >> I form of the speed distribution. 

Analysis: We shall consider first the stochastic model where the computations are simpler and essentially rigorous. 
The system now consists of N point particles in the unit 2D torus which move according to (jTJ) between collisions (with- 
out the term J 7 ,). In addition each particle independently has a (virtual) collision with a Poisson rate equal to A(qj)|v,| 
for some position dependent rate A(q) > 0, i.e., the weighted mean free path between collisions J Q A(q(t))|q(t)|di is 
an exponential random variable with mean one. The collision changes the angle which v makes with the x axis from 
9' to 9 according to some transition kernel K(9, 9') d9. The exact form of K will turn out not to matter as long 
as K(9,9') = K(6',9) and there is enough spreading to the direction of the velocity so that dqd9/2n is the unique 
invariant distribution for the system with one particle (N = 1) and E = 0. The scattering "closest" to that caused 
by collisions with fixed discs and the one we used in the simulations is the following: v changes to v' according to the 
rule 

v' = v - 2n(n • v) (3) 

where n is a unit vector in the direction of the momentum transfer from v to v'. The direction of n is chosen randomly 
with probability density — (n • v)/2, where v = v/|v|, subject to the constraint (n • v) < 0. 

The "master" equation describing the time evolution of the ./V-particle velocity distribution function is, for the 
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above rule, given by 

i— 1 z— 1 

AT 



i=l 



f A(<b)(vJ ' A ) [W(Q, V{, f; E) - W(Q, V, t; E)] dn 

./(v$-n)<0 



= AW + EBW + CW (4) 

where j = 3/U as in ©, 

Q = (qi,...,qjv) , V = (v 1 ,...,v i) ...,vjv) and V- = (vi, . . . , v<, . . . , v N ), (5) 

and ~v[ is given in terms of v, by (|3|). In the last term E is the magnitude of E, i.e., E = Ee for a unit vector e. 
Finally, C = an< ^ = SiLt ^ are the sum of collision and streaming terms which act independently and do 

not depend on E. We shall assume that A(q) is such that ^ has a unique steady state solution for every E ^ 0. A 
sufficient, but not necessary, condition is that C _1 < A(q) < C for some positive C. 

Let us consider now what happens when E is small. We note first that when E = the speed of any particle does 
not change with time but the collisions (deterministic or stochastic) randomize the direction of the velocity and the 
position of each particle. The distribution of speeds would then remain unchanged in time. When E is small the 
appropriate time scale for the change in the speed of the particle will be of order E~ 2 . Now on that time scale each 
particle will have undergone many collisions and so one may then assume that the direction of the velocity and the 
position of each particle will be uniformly distributed. We can thus expect to have an autonomous equation for the 
distribution of the speeds. Let us set v, = ?*j(cos(#j), sin(#j)) where r, = |vj| and the angle #j is taken with respect to 
the field direction which we can assume is in the a;-direction. Moreover we set R = (r\, . . . ,rjy) and = (6^, . . . , On). 
We then carry out a van Hove (weak coupling) limit [H, [l5|, rescaling of the time by letting t = t/E 2 and set 
W(Q, V, r; E) = W(Q, V, tE~ 2 ;E). W satisfies the rescaled equation 

cW(Q V,t;E) = E _ 2 + C ^ (Q Vj E) + E -l B w(Q, V, r; E) (6) 

OT 

We now assume that 

W(Q, V, r; E) = W {0) (Q, V, r) + EW {1) (Q, V, r) + E 2 W {2) (Q, V, r) + o(E 2 ) (7) 

This is a very strong assumption that can be better justified with a more detailed analysis Q . 
Substituting ([7]) into we get the following set of equations 

= (A + C)W^(Q,V,t) (8) 
= (^ + C)M? (1) (Q,V,r) +Bp)(Q,V,r) (9) 

aw(°)(Q,v,T) 



dr 



(A + C)W^(Q,V,t) + BW ( - 1 \Q,V,t) (10) 



Eq.© implies that W^(Q,V,t) depends only on R, i.e., (Q,V , r) = W(°)(R,t). Since BF(V) is orthogonal 
to the functions that depend only on R if F depends only on R, it follows that Sq^BW^ (R, r) = 0, where £q© 
is the average on Q and ©. Thank to our hypotheses on A and K we have that (A + C)F(Q, V) = if and only if 
F(Q, V) is a function of R alone so that W^(Q, V,r) = -(A + C)- 1 BW {0 \K, t) is well defined. We now insert 
this expression into (fTU)) and average over and Q. This does not effect the left hand side since depends only 
on R but it makes the first term on the right hand side vanish leading to 

^f' T) = -e & , Q B(A + crBwW(R,T). (11) 
which is indeed an autonomous equation for W^°'(R,t). 
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Equation can be written out explicitly as: 



^ dW^r) = gg^-[M -(R)^°)(R,r)] + £ ^(R)^°>(R,t)] (12) 



i=l j=l J i=l 

where the component of the N x N matrix M are given by 

m iM) = ± bt<mm _ 1 Sli _ n±a + ^ fv„ (13) 

fc=l ! fe=l 

and 

JV iV 

^( R ) = -0E^ + ^E^ 6* = ^"^ (14) 

fc=i fc=i 

The diffusion constant D in (fT2]) is just the integral of the velocity autocorrelation in the field direction e when the 
magnitude of the field E = 0, i.e., D = e • De, where 



1 f°° 1 f°° 

V=rn <Vl®Vl(t)> = T ^ T / ( Vl ® e ^+ Ci )* Vl ) A, 
l v i| Jo l v il Jo 



(15) 



and (•) is just averaging with respect to the uniform measure dqd8/ (2tt) that is stationary for E = 0. 

D is in fact the only term in (|12[) which depends on the collision kernel C in (j4]). For a spatially uniform and 
isotropic scattering, i.e. when A(q) is a constant and K(8',8) = K(9' — 8), we get D = DI with 

1 />27T />00 -i />00 />27T 

D = — dO [cos8 cos 8{t)]dt = — / (ft/ ddcos(8)e tCl cos(<9). (16) 
27r Jo Jo 2?r Jo Jo 

For the specific model used in (H|), D = 3/4. In the case of the deterministic billiards D will depend on the shape of 
the table. Note however that the NESS corresponding to the stationary solution of (|12[) is independent of D which 
really just sets a time scale (r ~ t/(DE 2 )). 
We note that 

M = SS* with Sy(R) = M?l (17) 



j 

which implies that (fP2"]) corresponds to a stochastic time evolution described by the Ito stochastic differential equation 

N 

dr t = -DA^R) dt + Y^ VDV2Sij(K)dBi, (18) 
i=i 

where Bi are N independent Brownian motions. One can in fact first derive (fT8"]) and then obtain (|12p [3l|. Using some 
general theory [§] this shows that the solution of ([12"]) is unique. 

Let now W(Q, V; E) be the stationary solution of (g]). Averaging W(Q, V; E) on Q and to get W(R] E) and then 
taking lim^^o VF(R; E) we get the stationary solution of ([12")) , W^°\ This is in fact the limit as r — > oo of W^(R, t) 
and is independent of D. To compute it we observe that if W(Q,V, t;E) solves (01 so does W(Q,V, i;E) = 
/i([/)W(Q, V, t; E) every positive function h. Moreover (Q| is invariant under the rescaling 

V^pV, t^p-H, E->p 2 E. (19) 

This suggests to look for Wo of the form 

W i0) (R) = h(U)F (R) (20) 

where Fo(pR) = p 2Ar_1 Fo(R) and /i(J7) assures that Wq has integral 1. With this assumption we get that Fo satisfies 
the equation 

ffI^ + !^-0 (21) 

^ \ n drf + Udn> 1 J 
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This equation can be easily solved and we get, when the initial state is such that U = N, 



W (0) (R) = -5{U-N) 



JV - 2 ™-i 

„3 



E^ 3 



i=i 



where Z is just the normalization 



Z 



N -, - 2 "- 1 AT 
3 



E 



iV 



lim f (r;N) = Cexp(-cr 3 ) 



where 



c = ^ 



2-r (i; 



0.2325 



= 0.5355 



Going beyond the limit fi^Owe find the first order correction (in E) to the stationary solution of (HJ): 

W(R, 0; E, N) = W {0) (R) + EW {1) (R, 0) + o(E), 

where 



JV 



W^(R, 0; N) = (^ + C)-W ( °)(R;iV) = J Fi(R)^r i c(q i) ^) 



(22) 



(23) 



To get the one particle marginal speed distribution fo(r; N) one has to integrate (|22|) over the variables r 2 , . . . , r^- 
We remark here that when TV — > oo we have 



(24) 



(25) 



(26) 



(27) 



with Fi(R) = (2iV-l) 



and 



c(qi,0i) 



e *(^+ c «) cos ^^= -(A+C,)" 1 ' 



(28) 



Ci and are the collision and streaming operators defined on the right hand side of (@|. Note that, for N — 1, the 
invariant solution to ((4]) is simply 



W(r, 6;E,1) = ^- + ^c(q, 9) + o{E) 
so that c(q, 0) is simply related to the N = 1 problem. 

Deterministic Billiard: The master (Liouville) equation for the deterministic model is given by 

dW d (Q,V,t) 



(29) 



at 



AW d + EBW d + C d W d 



(30) 



with C d W d representing the collisions with the fixed convex obstacles. In this case the stationary state is not absolutely 
continuous with respect to the Lebesgue measure, i.e. it will not have a smooth density Wd[Q, V;E), see @, 0- On 
the other hand it will have a smooth density W d (Q, V, t; E) for finite time t and we can consider the measure 



A4 )JV (dQ, dV) = W d (Q, V, t; E, N)dQ dV 



(31) 



For N = 1 we know that //^(dq, dv) = lim^oo /i^ 1 (dq, dv) exists and its projection on the energy surface |v| = 1, 
can be written as p,E,i{dq, dff) = dc\d0/ (2ir) + E5fj,i(dq, d0) + o(E). We note that the property p9[) remain true for 
any solution W d (Q, V,£;E) of (|30[) . The expansion can be generalized to the deterministic billiard by replacing 
c(q, 0) with <5//i(dq, d6). More precisely we have 



ME,/v(rfQ, <fV) = S(U - N) F (R)JQ dV + £F X (R) ^ rf <5/n(d qj , d6gdQ* dV 1 + o{E) 



(32) 
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FIG. 2: Comparison between the 1 particle speed marginals of the different tables with 5 particles. 

where dQ = Ili^li! with dq the normalized restriction of the Lebesgue measure to T 2 \obstacles, rfQ' = 
Ylj-ndqj, dV l = Ylj^dvj and -Fi(R) is defined after (f2"T)) . The above expression implies that /}at(g?Q, dV) = 
lmiE-j-o Ae.aK^Q, <~V) is absolutely continuous with respect to the Lebesgue measure on the energy sphere and de- 
pends only on the speeds. We cannot prove this statement rigorously but it is well verified by our numerical simulations 
involving the full billiard table or just a portion of it. 

Numerical results: We have concentrated our numerical simulation on the deterministic billiard system using the 
billiard tables depicted in Figure [TJ Table A is the same used in Q . Table B has the central obstacle moved down to 
break the symmetry but remaining rather close to Table A. Table C instead was chosen to be as asymmetric and as 
far from Table A as possible. 

We computed the one particle marginal of the speed distribution for all 3 tables. To obtain an accurate and reliable 
result we ran a very long trajectory recording the speed of particle 1 every time of the order of E~ 2 . In this way we 
can assume that the data we collected form a random sample from the distribution fo,d(r,N). This allows us to use 
the Kolmogorov-Smirnov test to check whether /o,d(?", N) = /o(r, N) described after (JUJ), see fiol H2I]. 

In figure[2]we plot the marginal distribution for all three table when E = 0.015625(cos(</>), sin(</>)) with <fr = ir/2 and 
N = 5 together with the theoretical prediction coming from (|22|) . In figure [3] we plot the similar results for N = 512 
and compare it with (|24| . The P-value of the KS test for the cases shown in these figures was greater than 23% giving 
a strong evidence that our hypothesis on the distribution of the observed data is indeed correct. A more extensive 
report on these simulations can be found at http://www.math.uab.edU/~khu/g/gt/speed.html 

To check whether the full distribution in (|32[) is valid for the deterministic billiard we chose N = 3 and fixed U = 3 
so that we can take r| = 3 — r\ — rf. Fixing S r = -\/3/30, i e t Aij(ri,r2) be the characteristic function of the square 
of side S r centered at ((i + 0.5)5 r , (j + 0.5)S r ). We ran 1000 trajectories of average length 2 • 10 s system time units 
for Table B with E = 0.04(cos(</>), sin(0)) with <fi = tt/3. Since it is very hard to collect enough statistic to properly 
resolve the full distribution, we did not attempt to decorrclate the the data nor to run any statistical test on the 
reliability of the following results. 

We used the above described trajectories to compute 

G c (i,j;k) = /i E ,3(A 4 j(Yi, r 2 ) cos(fc6>i)) G s (i,j;k) = #E,3(Ajj(ri, r 2 ) sin(fc6>i)) 

and 

L c (i,j;k) = AE,3(Aij(ri,r 2 )cos(fc6)i)xA(qi)) L s (i,j;k) = £E,3(Ai,j(ri, r 2 ) sin(k6i)xA(<li)) 
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FIG. 3: Comparison between the 1 particle speed marginals of the different tables with 512 particle. 



with k = 1,2,4,8. Here Xa(<i) is the characteristic function of the lower quarter of the table, i.e. Xa(q.) = 1 if 
q x < 0.5 and q y < 0.5 and otherwise. Finally the we ran 200 trajectories of length almost 10 9 system unit time for 
an identical system with N = 1 to compute g c (k) = <5/ii(cos(fc#)) and l c (k) = 5/xi(cos(fc#)xA(q)) and the relative sine 
terms. In this way we can check (|3"2"j) with no parameters to be fitted. Observe moreover that from (p?2"j) we have 

G c (i,j; k) = a c (k)L c (i,j; k) 

with a c {k) = g c (k)/l c (k), and similarly for G s . In figure 2 we plot G s (i, 20; 1) with a s (l)L(i, 20; 1) and the theoretical 
prediction. We have selected to plot only j = 20, corresponding to T2 = 1.184, for better readability. On the horizontal 
axis we have T\ instead of i for better readability. As one can see, the agreement is very good. Similar agreement 



is found for the other k. The interested reader may visit http : / / www . math . uab . edu/~khu/ g/ gt/ speed . html| where 
more plots of the results of these simulations can be found. 
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